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ABSTRACT 


A dominant height growth model and a site index model were developed for Calabrian pine (Pinus brutia 
Ten.) in central Cyprus. Data from 64 stem analysis in 32 temporary plots, where Calabrian pine was the 
only tree species, were used for modeling. The plots were selected randomly in proportion to two site 
types. Four difference equations were tested. The evaluation criteria included qualitative and quantitative 
examinations and a testing with split data. The difference equation of Korf showed the best results for 
all data. An analysis of the height growth patterns among sites - as these were defined from the selected 
equation - was made in order to study the behavior of different site index curves. Results indicated the 
validity ofa common height growth model for the two sites. In spite ofthe irregular height growth pattern 
observed in Calabrian pine, the model obtained allows us to classify and compare correctly Calabrian pine 
stands growing at different sites. 
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1. Introduction 


The species Pinus brutia Ten. (Calabrian pine) is a characteristic 
species of the eastern Mediterranean [12,26,27,34,64,83]. P. brutia 
is widely extended in Turkey and far Eastern Greece, secondarily in 
the Crimea, Caucasus coast, Azerbaijan, Iraq, Syria, Lebanon, Crete 
and Cyprus [54]. 

It is an important species for multi-purpose forestry; construc- 
tion, industry, carpentry, firewood and pulp [7,10,61 ]. Additionally, 
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its post-fire regeneration ability makes it a unique forest species 
within the fragile Mediterranean ecosystems [68]. 

P. brutia generally occurs in the Mediterranean-type climate of 
hot and dry summers and mild and rainy winters. This species 
occurs most abundantly in the semi-arid and sub-humid zones 
[64]. The species has strict rainfall requirements, being absent from 
the arid bioclimatic zone and rare in the lower semi-arid zone, 
but widespread in the humid zone [63]. We find Calabrian pine in 
zones with a mean annual precipitation between 400 and 2000 mm 
[36,45]. P. brutia is well adapted to the Mediterranean-type cli- 
mate in several physiological and morphological characteristics 
and is a drought resistant species [17,19,33,44,45,53]. The adapt- 
abilities of its provenances to drought vary [32], but P. brutia 
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generally achieves its optimum growth in rainy regions of more 
than 900-1000 mm mean annual precipitation. Its typical elevation 
range is between 0 and 1500 m above sea level [45]. The geograph- 
ical isolation among populations implies the existence of several 
races, with different ecological characteristics and behavior [54]. 

Calabrian pine is a fast growing tree species; its mean annual 
increment in plantations can be over 10 m?/ha (at a spacing of 
4.5 m?; 2222 trees per hectare, site class I, site index 20 m) between 
the ages of 25 and 35 years [35]. New stands of P. brutia generally 
begin after fires, although it can also naturally regenerate without 
fire [26,27,46,47,48,76,77]. Several adaptations greatly contribute 
to its post-fire regeneration [20,46]. Seedlings develop rapidly 
growing tap roots while the stems grow comparatively slowly in 
height, so that seedlings can have 65 cm tap roots within five to six 
months after germination [26,27,46]. Early stem growth is greater 
at higher elevations, while average root growth is the opposite [77 ]. 
The early rapid growth of P. brutia continues during the following 
years; the mean annual increment can be over 10 m?/ha in appro- 
priately spaced plantations between the ages of 25 and 35 years, as 
described earlier (site class I, site index 20 m [35]. 

Calabrian pine is the main forest species in Cyprus growing from 
sea level up to the 1600 m altitude [79]. Its forests represent about 
90% of the Island’s forest area, which totals at 173,182 ha [5]. Most of 
theseforests were classified currently unexploitable after a recently 
carried out inventory while 41,399 ha were classified as currently 
exploitable [5]. Fellings are confined in the latter. The annual wood 
output for the periods 1988-1992, 1998-2002 and 2003-2007 was 
56,100 m3, 27,633 m? and 14,384 m3, respectively [5]. This wood 
is utilized locally for the production of sawn timber, chipboard, 
box-shooks, fuelwood, charcoal and other minor wood products 
[79]. Sixteen thousand eight hundred thirty five hectares of the P. 
brutia forests have been selected and preliminarily classified either 
as National Forests Parks to provide amenities to the public, or as 
Nature Reserves for nature conservation [5]. 

The past management (silvicultural system) of pine forests in 
Cyprus proved to be unsuitable, resulting in serious ecological, 
economic and social constraints [79]. This may put in danger 
the survival of these stands in the long term [43]. The degrada- 
tion of pine forests in Cyprus has resulted in drastic reduction of 
annual wood output from 56,100 m? (for the period 1988-1992) to 
14,384 m3 (for the period 2003-2007) [5]. The recognition of this 
situation and the increasing interest in using these stands for either 
direct production, or indirect production (environmental preserva- 
tion) justifies the need for a sustainable management of Calabrian 
pine stands. Considering the high environmental and silvicultural 
variability of P. brutia stands [43], it is necessary to typify and char- 
acterize them in order to optimize their management. 

Estimating forest productivity is both necessary for effective for- 
est management and useful for evaluating basic site conditions for 
ecological field studies. Site quality is therefore influenced by fac- 
tors such as available light, heat, moisture, and nutrients, along 
with other soil characteristics such as soil depth and aeration [29]. 
Although it would be best to directly measure and predict these 
factors, they require precise measurements that may be difficult 
to extrapolate across scales. Therefore, indirect methods for evalu- 
ating site quality are more frequently used in forest management 
[56,67,72]. 

Site index, defined as dominant height at some fixed base age, 
is one of the most commonly used indicators of site productivity 
because there is a close correlation between volume and dominant 
height growth [72]. Many mathematical functions are available 
to model dominant height growth. Desirable characteristics for 
growth functions are: polymorphism, existence of inflection point 
and horizontal asymptote, logical behavior, right theoretical basis, 
base-age invariance, and parsimony [9,14,65]. These requirements 
depend on both the construction method and the mathematical 


function used to develop the curves [65]. Among other meth- 
ods for site index curve construction [38], the algebraic difference 
approach presents certain advantages: short observation periods 
can be effectively used, and the structure of equations is base-age 
invariant [9,14,15,42]. 

Modeling dominant height growth for Pinus sp. in the Mediter- 
ranean region has received some attention. Río and Montero [58] 
developed site index curves for Pinus sylvestris in Spain, Calama 
et al. [67] for Pinus pinea in Spain, Bravo-Oviedo et al. [2] for Pinus 
pinaster in Spain, Hatzistathis et al. [4] and Kitikidou et al. [41] for 
P. brutia in Greece. 

Shater et al. [85] developed dominant height, individual-tree 
diameter increment, tree height, and self-thinning site curve mod- 
els for P. brutia stands in Syria. In this work, the Chapman-Richards 
line was selected to show the development of dominant height 
on average site. Palahí et al. [57] developed and compared two 
sets of models (site index models for even-aged and uneven-aged 
stands) which enable tree-level simulation of the development 
of pure and mixed stands of P. brutia in North-eastern Greece. 
Adamopoulos et al. [73] examined the impact of site qualities on 
P. brutia wood characteristics in North-eastern Greece, while Kav- 
gaci et al. [6] studied the post-fire long-term regeneration of the 
species in South-western Turkey through a floristic survey. Aert- 
sen et al. [80] compared and evaluated five different modeling 
techniques (multiple linear regression, classification and regres- 
sion trees, boosted regression trees, generalized additive models 
and artificial neural networks) to model site index of homogeneous 
stands of important tree species of Turkey (P. brutia among others) 
by using environmental variables (soil, vegetation and topograph- 
ical variables) as predictors. In another study, Aertsen et al. [81] 
compared the performance of the same five modeling techniques 
forthe prediction of forest site index (P. brutia stands among others) 
in two contrasting ecoregions (lowlands in Belgium and Mediter- 
ranean mountains in Turkey). Artificial neural network models 
have been developed to estimate volume of dominant Calabrian 
pine trees [49,51]. Finally, Diamantopoulou et al. [50] applied 
Principal Components Analysis and Hierarchical Cluster Analysis 
in order to interpret the behavior of several quantitative vari- 
ables of reforestation Calabrian pine trees (dominant height among 
others). 

The main goal of this study was to develop a dominant height 
growth model for P. brutia growing in two different site types in 
central Cyprus, from which a site quality model would be created. 
To achieve this objective, the variability in dominant height growth 
patterns and the differences among the sites were analyzed. 


2. Materials and methods 
2.1. Data set 


The study was conducted in central Cyprus (34°56N, 33°17E), 
in middle elevation (300-750 m a.s.l.) forest areas. This area cov- 
ers about 14,000ha. The annual rainfall and the mean yearly 
temperature (for 1991-2000), according to the closest meteoro- 
logical station (Kornos, 34°55N, 33°24E, 370 m a.s.1.), are 477.1 mm 
and 18.4°C respectively. The geological substrate of the study 
area, belongs to volcanic sequence diabase dykes with pillow lava 
screens [30]. The soil texture is mainly sandy loam to sandy clay 
loam (Petrou and Milios [62], data from Forest Department of 
Cyprus). 

The study area was formerly consisting mainly of vine and olive 
yards. The majority of those fields were abandoned during the last 
century and nowadays they have been recolonized mainly by P. 
brutia trees which is the dominant species in the area. A few olive 
yards are still cultivated. Moreover three built up areas are located 
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Table 1 
Summary statistics of sampled trees per site. 
Site No. of trees Mean T SDT Tmin T max Mean H SDH H min H max 
] 32 27.69 4.12 24 37 13.06 2.04 9.90 18.80 
2 32 90.09 743 76 108 8.93 1.67 5.70 12.30 


H (m)is total tree height at age T (years). 


in the study area. These are the villages of Lythrodontas, Mathiatis 
and Kapedes. 

The P. brutia trees form stands with different stand structures 
and variable densities in the wider study area. In some cases the 
P. brutia trees create small single cohort stands [18]. The P. brutia 
trees in those stands have a ground cover percentage of at least 75%, 
forming dense stands. Also the P. brutia trees can be classified to 
dominant, co-dominant and intermediate according to their crown 
class. [18,22]. 

Particularly the abovementioned single cohort stands are 
located in two different site types as regards their soil depth (site 
types A and B). The soil depth in the plots of site type A is over 
100 cm, while the corresponding depth in the plots of site type B is 
about 30 cm. The soil depth was measured at 2 randomly selected 
points in each plot. Thirty-two plots of 100 m? (10 x 10m) were 
randomly established in both site types (16 in site type A and 16 in 
site type B). From each plot one dominant and one co-dominant tree 
was felled, according to the stratified random sampling method, 
for stem analysis purposes, resulting in 64 trees (32 in each site 
type). From each tree, cross-sectional discs were cut and removed 
from the stump height («0 m), at breast height (1.3 m) and every 
1m, up to the bole. The last disc was collected from the 3 cm-bole 
diameter. In each cross-sectional disc the number of annual growth 
rings was counted using the LINTAB system of RinnTech and the 
program TSAP-Win [25]. In stem analysis for the height calcula- 
tion the modified version of Carmeans' algorithm was used [40,82]. 
The data set was composed of 3769 recordings from the 64 sam- 
pled trees. The data summary of the sampled trees is presented in 
Table 1. 

Data were split in two groups. The first group (80% of the trees 
- 2905 recordings from 51 trees) was used as fitting data, and the 
second group (20% of the trees - 864 recordings from 13 trees) was 
used as validation data. This method for testing the performance of 
a proposed model is known as cross-validation [21,29,52,59,60,70]. 


2.2. Candidate functions 


An algebraic difference approach has been used, since it shows 
better properties and performance than analogous fixed-base-age 
equations [13,15]. Four models were selected from those most com- 
monly used in forest research as candidate functions to model 
dominant height growth (Table 2). The first two polymorphic func- 
tions derive from the Bailey-Clutter function [66], the next two 


Table 2 
Candidate models for dominant height modeling. 


from the Korf function [78], the next is the McDill-Amateis func- 
tion [55], and the final two polymorphic functions derive from the 
Chapman-Richards function [24]. 


2.3. Data structure and model fitting 


To fit an algebraic difference equation expressed in the general 
form of H5 =f (H4, Ti, T2), different data structures defined in Bor- 
ders et al. [8] can be used. These data structures are relevant in any 
increment-based modeling. The data chosen for fitting the differ- 
ent functions included all the possible combinations of height-total 
age pairs for a tree (all possible growth intervals) [8,28,39,67,74]. 

The fittings were carried out using the nonlinear procedure on 
the SPSS software [71]. The Levenberg-Marquardt iterative method 
was selected because it is the most useful when the parameter 
estimates are highly correlated [84]. 

The autocorrelation derived from using stem analysis data 
was barred by applying the Goelz and Burk [39] correction. First, 
each function is fitted following ordinary non-linear least squares 
regression and the error term ej, residual from estimating H; using 
Hj, is expanded following an autoregressive process: 


eij = p&i-1,j + YEij-1 t €ij 


where: o-autocorrelation between the current residual and the 
residual from estimating H; , using Hj as a predictor variable; 
y 7 relationship between the current residual and the residual form 
estimating H; using Hj. , as a predictor variable; ej ; -independent 
errors with mean zero and constant variance 7”. The model param- 
eters are then obtained by fitting the expanded function. The 
autocorrelation parameters vary the weight of each observation by 
reducing the residual proportional to a previous residual. Neither 
of the autoregressive parameters p nor y are used for field appli- 
cations of equations because the errors e£; , ; and £; j 1 cannot be 
observed without stem analysis [72]. 


2.4. Model selection criteria 


Athree-step procedure was used to evaluate and select the most 
appropriate model, which included qualitative as well as quantita- 
tive examinations. The first step was to evaluate the model fitting 
statistics based on seven model performance evaluation criteria 


Model Name Height-age equation Free parameter Difference form 

1 Bailey-Clutter [66] H-a,*aj;T6 a Hy = He 02 Tp) (1.1) 
a H5 = e^i «(n H4-aiXT2 /T; 3 (1.2) 

2 Korf [78] H =ayel%2/T3) a Hy = ai(Hi /a; "1/2" (2.1) 
a3 Aa = a Tay ain; Tin (a, /H, J9 72 (2.2) 

3 McDill-Amateis [55] H- THn/TS) az EEG T (3) 

4 Champan-Richards [24] H -aj[1 — e20" ] a H =H, (1R)^ (4.1) 
a Hy =a, f1- ji- (8) Pe (42) 


Symbols: H is total tree height at age T; H and H; are heights at age T? and T;, respectively; a; are model parameters and e is the base of the natural logarithm In. 
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Table 3 
Model performance evaluation criteria (fitting and testing procedures). 


Criterion Symbol Formula Optimum value 
n 
Mean residual MRes 1 BH 0 
il 
Rol 
, Uli -H) 
Variance ratio VR DE 5 1 
pL 
n 
Residual mean of squares RMS 1 DD 0 
ii 
n " 
S SR] 
Absolute mean residual AMRes = 0 
i=1 
Liam 
Coefficient of determination-model efficiency R?-MEf 1 i = 1 
Linear regression a, B Hi =a + BA; + £i a=0, B=1 


Symbols: fl; is the ith estimated value; H; is the ith observed value; n is the number 


described by Amaro et al. [1] (Table 3), and select those equations 
who appeared to be the best. 

In step two, the characterization of the model error was ana- 
lyzed, based on the independent data set testing defined in Section 
2.1 (864 recordings from 13 trees). For validation purposes, evalu- 
ation criteria applied in the first phase were also calculated. 

Finally, the correctness of the theoretical and biological aspects 
of the models was assessed in step 3. 


2.5. Comparison of height growth models among sites 


Once the best function had been selected, the differences in the 
dominant height growth models for the different sites were com- 
pared using both the full and the reduced models. The reduced 
model corresponds to completely different sets of parameters for 
different sites - as these are defined from the selected function - 
while the full model corresponds to the same set of parameters for 
all the sites combined. 

Two tests for detecting simultaneous homogeneity among 
parameters were used: the Bates and Watts non-linear extra sum 
of squares F test [74,75] and the test proposed by Lakkis and Jones, 
in Khattree and Naik [69], to compare the differences in site index 
models between sites. These tests are frequently applied to analyze 
differences among different regions [37,67,75]. 

Besides the full and reduced models, the sum of squares error 
(SS) is necessary to calculate both tests. This error was calculated 
as follows: 


where: n=number of observations for each tree; m=total number 
of trees. 
The F-test is applied using the following equation: 


(SS, — SS¢)/(dfr — dfr) 
(SSp)/(dff) 


where: SS; and SS, = error sum of squares for full and reduced model 
respectively; df and df; - degrees of freedom for full and reduced 
model respectively. F follows an F-distribution. 


of observations; p is the number of model parameters. 


Table 4 
Fit statistics for each model. 


Model MRes VR RMS AMRes MEf Linear regression 
a B 

M1 4.0092 0.0632 25.5877 4.0092 -1.0147 2.0971 2.4969 

M2 —0.0145 0.3855 7.6554 0.0145 0.3972 .-0.0954 1.0153 

M3 —0.0359 0.3611 7.7283 0.0359 0.3915 .-0.2602 1.0421 


The L statistic used in the Lakkis—Jones test is defined as: 


m/2 
a eel (9) 
SSr 


where: SS; and SS, = error sum of squares for full and reduced model 
respectively; and m = total number of trees. If homogeneity exists 
among the model vectors of parameters f, the distribution of the 
statistic —2In(L) converges in probability to a Pearson x? distribu- 
tion, with v degrees of freedom, where v is equal to the difference 
between the number of parameters estimated in the full model and 
the reduced models. 


3. Results 


Performance results of all models are shown in Table 4. All the 
parameter estimates for all the functions were significant at an 
o level of 5%, except for those of the Champan-Richards function 
(Model 4) which did not meet the convergence criterion. The mean- 
ing of the negative value of MEf for the Model 1 is that, the mean 
value of observed heights gives a better estimate than the expected 
values resulting from this model [1]. The Korf [78] model resulted 
to be the best model. Differences between M2 and M3 were very 
small, although the analysis of the fit statistics revealed that the 
Korf function generally resulted in slightly lower values for MRes, 
RMS and AMRes as well as higher efficiency. The shape of the curve 
of the Model 2 is independent of the autocorrelation correction 
(Fig. 1 shows both fits of M2), showing practically the same curve. 
Regression coefficients for M2 are calculated as: 


T0.648) 


H = 10.924e(-5:876)/( 


H = 10.9266 (-5:875)/(19947) 


without and with autocorrelation correction, respectively. 
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Table 5 
Testing statistics for each model. 
Model MRes VR RMS AMRes MEf Linear regression 
a B 
M1 3.5235 0.0892 19.8586 3.5235 —0.8877 2.2185 2.1629 
M2 —0.0015 0.4404 5.8193 0.0015 0.4468 —0.0353 1.0073 
M3 —0.0092 0.4393 5.8113 0.0092 0.4476 —0.0530 1.0094 
94 30 
84 
25 4 
74 
E «6! E 20 4 
E $ 
$5 3 
E E as] 
È 4] E 
E 5 
&* a 10] 
—— Model 2 
24 * - - Model 2 with autocorrelation correction 
54 
14 MT ——— Model 2 - Site | 
m cu - ++ -Model 2 - Site II 
0 at 
0 4 8 12 16 20 24 28 32 36 40 44 48 0 10 20 30 40 50 60 70 80 90 100 110 


Total age (years) 


Fig. 1. Shapes of the curves of Model 2 resulting from both autocorrelation correc- 
tion and not autocorrelation correction fits. 


The functions M2 and M3 were again those which per- 
formed better with the testing data (Table 5). Once more, 
Chapman-Richards function did not converged, and the Models 2 
and 3 performed similarly. 

The shapes of the curves of Model 2 seem to be biologically 
reasonable, which prevent unrealistic height predictions when 
extrapolating the function beyond the range of the original data. 
Model 2 was selected, based in the good results in the performance 
criteria defined in Table 3 with both fitting and testing data. 

Assuming the suitability of the Model 2, it would be interest- 
ing to analyze the dominant height growth pattern among sites. 
Figs. 2 and 3 represent the Model 2 adjusted for each site, as these 
are defined from Model 2. Site I shows a very rapid dominant height 
growth pattern, especially for younger ages. 

The non-linear extra sum of squares F test and the Lakkis—Jones 
test revealed that, the null hypothesis of parameter homogeneity 
was acceptable in the full model for all sites. The reduced model 
for site I did not converged (Table 6). On the other hand, the fit 
statistics obtained with full and reduced models, in all sites and 


Total age (years) 


Fig. 2. Site index curves for Pinus brutia in central Cyprus, using the Korf function 
(Model 2). 


sites I and II (productive and less productive) separately, were very 
dissimilar (Table 7). According to these results, the total model (all 


sites together) can be selected. In our attempt to apply the selected 


M2 


- Model 2 - Site | 


Height (m) 


Age (years) 


Fig. 3. Site index curves for Pinus brutia in central Cyprus using the Korf function 
together with the stem analysis. 


Table 6 
Land F statistics. 
M2 df full model SS; df reduced model SS, Parameter L Parameter F 
All sites 2902 22216.522 2902 22216.522 1.000 0.000 
Site I Did not converged 
Site II 744 1461.900 2155 3337.704 0.0000* 0.6766 
Significant L-values and F-values are marked with *. 
Table 7 
Fit statistics of the full and reduced models for the two sites. 
M2 N Model MRes VR RMS AMRes MEf Linear regression 
o B 
All sites 2905 Full —0.0145 0.3855 7.6554 0.0145 0.3972 —0.0954 1.0153 
Reduced —0.0142 0.3855 7.6560 0.0142 0.3972 —0.0943 1.0151 
Site I 747 Full 0.0333 0.9271 1.9649 0.0333 0.9035 0.1135 0.9873 
Reduced Did not converged 
Site II 2158 Full 4.8979 0.0514 36.7628 4.8979 —0.8063 0.3312 4.1685 
Reduced 1.2765 1.2269 23.3246 1.2765 —1.4401 5.7317 —0.2036 


1328 K. Kitikidou et al. / Renewable and Sustainable Energy Reviews 16 (2012) 1323-1329 


model to sites I and II separately, the reduced model for site I did 
not converged, while both full and reduced model for site II gave 
results far beyond ideal. 


4. Discussion and conclusion 


Accurate estimates of forest productivity are needed for sustain- 
able forest management in order to determine annual allowable cut 
and rotation period, and to make tree species selection decisions 
[11]. When the time needed to produce economically feasible and 
ecologically sustainable crops is unknown and existing growth and 
yield information is limited, developing site index models helps to 
investigate growth performance [11]. 

This study presents a site index model for the P. brutia stands 
in central Cyprus. The algebraic difference approach used to derive 
dynamic height functions of dominant trees assures a high degree 
of robustness in applications [15]. According to fit statistics, the 
growth function of Korf [78] was chosen to explain the height 
growth pattern of this species. The statistic criteria with testing 
data showed that the growth function of Korf was still the best pat- 
tern. Testing using separate intervals (sites) produced worse results 
that working with all possible intervals (sites), so it seems advisable 
to analyze all data together, especially in the testing phase, in view 
of the fact that the use of separate intervals is not very common in 
site index estimation [23]. 

Comparing the results of the selected model (M2) among sites, 
there were some differences in error measurements and model 
efficiency values. Similarly, parameter estimates varied consider- 
ably among sites. This could be due to little data available for an 
adequate estimation [16] as a result of the absence of old stands, 
probably because of their previous land use as vine and olive yards. 

An effort was made to ensure that sampling occurred through- 
out the entire range of Calabrian pine stands in the study area. 
Within this range, it seems that an intermediate site is missing 
(Fig. 3). Although Calabrian pine is found in intermediate sites in 
the study area, the pine trees are not forming so dense stands. In 
this sense, the general population of pine stands was sampled, with 
canopy cover at least 7525 of the ground, as mentioned in Section 
2.1. Thus, the model that was developed provides site index esti- 
mates that will be applicable throughout most of the range of rather 
dense pine stands in the study area. 

In spite of the apparent differences observed also in the graph- 
ical comparison of the site index curves (Figs. 2 and 3) obtained by 
fitting the selected function to each site, the statistical tests did not 
reject the null hypothesis of equality of height growth patterns. 

The application of a single height growth model for all the stud- 
ied area, can help us classify the quality of Calabrian pine stands and 
decide what measures can be taken regarding the different silvicul- 
tural treatments [37]. Nowadays, due to the state of degradation in 
which many of P. brutia stands in Cyprus have fallen [79], itis neces- 
sary to establish priority areas where resources should be invested 
and this site index model may provide a key tool for this purpose. 
A common dominant height growth model could also simplify the 
development of a site index model based on ecological variables 
[3]. Although it is sometimes necessary to stratify the study area 
when looking for relationships between site index and ecological 
variables, the analysis is facilitated with a single dominant height 
growth pattern [29,31]. 

As regards comparison of P. brutia growth patterns in the 
Mediterranean region, if the obtained curves in this work are com- 
pared with the site index model for Greek P. brutia stands in Thasos 
island [41], stands in central Cyprus show higher growth rates at 
younger ages and a larger gap between the best and the worst site 
qualities. This could be due to different soil depth of the two areas 
or due to different genetic background. Additionally, the large sam- 
pled area of this study and the good results obtained in different 


sites (fitting and testing data) makes this site index model a good 
option for classifying site qualities of P. brutia stands in Cyprus. 
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